library(lmtest)
library(car)
library(ebal)

## set this to the appropriate directory
setwd("C:\\Papers\\SF_natural_experiment\\Replication_Files")
options(scipen = 999) 

SF.data = read.csv("SF_ne_data.csv", header=TRUE)

eb.out = ebalance(SF.data$new.pp, subset(SF.data, select=c(quake_disruption, pct_over65, pct_college, pct_black, pct_asian, pct_female, med_income, consolidated)))
SF.data$weights = 1
SF.data$weights[SF.data$new.pp==0] = eb.out$w

## entropy weighted model ##
reg.w = lm(turnout.diff ~ new.pp + quake_disruption + pct_over65 + pct_college + pct_black + pct_asian + pct_female + med_income + consolidated, data=SF.data, weights=weights)
summary(reg.w)
bptest(reg.w)
cov = hccm(reg.w, type="hc1") 
weighted.results = round(cbind(reg.w$coefficients, sqrt(diag(cov)), reg.w$coefficients - 1.96*sqrt(diag(cov)), reg.w$coefficients + 1.96*sqrt(diag(cov))), digits=3)
colnames(weighted.results) <- c("Coeff.", "SE", "2.5th Pct.", "97.5th Pct.")
print(weighted.results)
write.csv(weighted.results, "weighted_results.csv")

## placebo test using adjacent precincts and excluding treated ##

SF.data.placebo <- subset(SF.data, new.pp==0)
eb.out.p = ebalance(SF.data.placebo$placebo, subset(SF.data.placebo, select=c(quake_disruption, pct_over65, pct_college, pct_black, pct_asian, pct_female, med_income, consolidated)))
SF.data.placebo$weights.p = 1
SF.data.placebo$weights.p[SF.data.placebo$placebo==0] = eb.out.p$w

## entropy weighted model
reg.w.p = lm(turnout.diff ~ placebo + quake_disruption + pct_over65 + pct_college + pct_black + pct_asian + pct_female + med_income + consolidated, data=SF.data.placebo, weights=weights.p)
summary(reg.w.p)
bptest(reg.w.p)
cov.p = hccm(reg.w.p, type="hc1") 
weighted.results.placebo = round(cbind(reg.w.p$coefficients, sqrt(diag(cov.p)), reg.w.p$coefficients - 1.96*sqrt(diag(cov.p)), reg.w.p$coefficients + 1.96*sqrt(diag(cov.p))), digits=3)
colnames(weighted.results.placebo) <- c("Coeff.", "SE", "2.5th Pct.", "97.5th Pct.")
print(weighted.results.placebo)
write.csv(weighted.results.placebo, "weighted_results_placebo.csv")


## unweighted models -- not reported (but useful to see results are similar to entropy weighted model) ##
reg.uw = lm(turnout.diff ~ new.pp + quake_disruption + pct_over65 + pct_college + pct_black + pct_asian + pct_female + med_income + consolidated, data=SF.data)
summary(reg.uw)
bptest(reg.uw)
cov.uw = hccm(reg.uw, type="hc1") 
unweighted.results = round(cbind(reg.uw$coefficients, sqrt(diag(cov.uw)), reg.uw$coefficients - 1.96*sqrt(diag(cov.uw)), reg.uw$coefficients + 1.96*sqrt(diag(cov.uw))), digits=3)
colnames(unweighted.results) <- c("Coeff.", "SE", "2.5th Pct.", "97.5th Pct.")
print(unweighted.results)
write.csv(unweighted.results, "unweighted_results.csv")

## unweighted placebo test using adjacent precincts
reg.uw.p = lm(turnout.diff ~ placebo + quake_disruption + pct_over65 + pct_college + pct_black + pct_asian + pct_female + med_income + consolidated, data=SF.data.placebo)
summary(reg.uw.p)
bptest(reg.uw.p)
cov.uw.p = hccm(reg.uw.p, type="hc1") 
unweighted.results.placebo = round(cbind(reg.uw.p$coefficients, sqrt(diag(cov.uw.p)), reg.uw.p$coefficients - 1.96*sqrt(diag(cov.uw.p)), reg.uw.p$coefficients + 1.96*sqrt(diag(cov.uw.p))), digits=3)
colnames(unweighted.results.placebo) <- c("Coeff.", "SE", "2.5th Pct.", "97.5th Pct.")
print(unweighted.results.placebo)
write.csv(unweighted.results.placebo, "unweighted_results_placebo.csv")

